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Abstract. Modelling gravity is a fundamental problem 
that must be tackled in iV-body simulations of stellar 
systems, and satisfactory solutions require a deep under- 
standing of the dynamical effects of softening. In a previ- 
ous paper (Romeo 1997), we have devised a method for 
exploring such effects, and we have focused on two ap- 
plications that reveal the dynamical differences between 
the most representative types of softened gravity. In the 
present paper we show that our method can be applied 
in another, more fruitful, way: for developing new ideas 
about softening. Indeed, it opens a direct route to the 
discovery of optimal types of softened gravity for given 
dynamical requirements, and thus to the accomplishment 
of a physically consistent modelling of disc galaxies, even 
in the presence of a cold interstellar gaseous component 
and in situations that demand anisotropic resolution. 

Key words: gravitation - methods: analytical - methods: 
numerical - galaxies: evolution - galaxies: kinematics and 
dynamics - galaxies: spiral 



1. Introduction 

A^-body simulations of disc galaxies rely on the use of 
softening. This artifice removes the short-range singular- 
ity of the gravitational interaction, which is dynamically 
unimportant and computationally troublesome, whereas it 
leaves the long-range behaviour of gravity unchanged. But 
softening is also a critical factor in simulations. It controls 
their quality and can affect their result on scales much 
larger than the softening length. Its dynamical effects are 
further exacerbated in the presence of a cold interstel- 
lar gaseous component and in situations that demand 
anisotropic resolution. Thus softening poses a dynami- 
cal problem of special concern, which should be probed 
carefully and in detail (e.g., Hernquist & Barnes 1990; 



Pfenniger & Friedli 1993; Romeo 1994, hereafter Paper I; 
Romeo 1997, hereafter Paper Iltl; and references therein). 

In Paper I, we have investigated how faithful simula- 
tions are. In particular, we have concluded that the stan- 
dard way of introducing softening in the presence of stars 
and cold interstellar gas is definitely unsatisfactory in sev- 
eral regimes of astrophysical interest. It is so because im- 
portant small-scale instabilities of the gaseous component, 
e.g. those peculiar to star-formation processes, are sup- 
pressed just as unphysical noise of the stellar component. 
Faithfulness requires an appropriate introduction of two 
softening lengths, one for each component, and also a rig- 
orous specification of the star-gas gravitational interac- 
tion. 

In Paper II, we have devised a method for exploring the 
dynamical effects of softening. As a major result, we have 
shown how to choose the softening length for optimizing 
the faithfulness of simulations to the Newtonian dynam- 
ics. Then we have focused on two applications that reveal 
the dynamical differences between the most representative 
types of softened gravity. In particular, we have concluded 
that it is desirable to improve the current way of introduc- 
ing anisotropic softening. We need a clearer decoupling of 
the resolution parallel and perpendicular to the plane, and 
also more natural planar and vertical softening lengths. 

In the present paper, which completes our planned re- 
search work about softening, we propose an innovative 
solution to the problem. The understanding of galactic 
and extragalactic astrophysics is at a crucial stage. Un- 
solved problems are viewed in new perspectives, which 
promise major revisions of knowledge (see, e.g., Blitz & 
Teuben 1996; Block & Greenberg 1996). Recent investi- 
gations suggest, for instance, a more enigmatic interplay 
between stellar disc and bulge/halo (e.g., Lequeux et al. 
1995), a clearer relation between cold gas and dark mat- 
ter in spiral galaxies (e.g., Pfenniger et al. 1994; Pfen- 
niger & Combes 1994; Combes & Pfenniger 1997), and a 



Sections and equations of that paper are denoted by the 
prefix II. 
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closer connection between the fractal structures of the in- 
terstellar medium and of the universe (e.g., de Vega et al. 
1996, 1998). The implications are clear: modelling grav- 
ity in iV-body simulations of disc galaxies should offer a 
flexible interface with such a progress. Our solution is to 
optimize the fidelity of simulations to given dynamical re- 
quirements. How do we apply this idea in practice? 

1 . We impose the requirements in the wavenumber space 
since this is the natural dynamical domain of gravity, 
as Pfenniger & Fricdli (1993) have previously empha- 
sized. 

2. We identify the softening length with the characteristic 
dynamical scale length. 

3. Then we invert part of the method of Paper II, and 
the result is the optimal type of softened gravity that 
satisfies those dynamical requirements. 

Our application covers both 2-D and 3-D modelling. 
The basic cases arc extended to more complex situations 
through recipes for implementing star-gas and anisotropic 
softening, which have already been motivated (cf. discus- 
sions of Papers I and II). Last but not least, each descrip- 
tion is complemented by an example that leaves room for 
creativity. 

The present paper is organized as follows. The appli- 
cation is shown in Sects. 2 and 3 (see also Appendix A), 
and proceeds as in the previous discussion. Comments on 
related works concerning softening arc made in Sect. 4. 
The conclusions and perspectives are drawn in Sect. 5, 
where we present our three papers about softening in a 
more unified view and emphasize their potentially strong 
impact on galactic dynamics. 



2. 2-D modelling 

2.1. Inverting part of the method of Paper II 

The method of Paper II allows determining the dynamical 
response of the model to a given type of softened gravity. 
The basic quantity that describes this response is the re- 
duction factor iS(|fc|s) defined in Eq. (II-4), where k has 
the meaning of a radial wavenumber and s is the softening 
length. In a few words, this is the factor by which soften- 
ing reduces the dynamical contribution of self-gravity. The 
behaviour of 5(|fc|s) provides the following information: 

• At small \k\s, it shows how significantly the stabil- 
ity and collective relaxation properties are affected on 
large scales. Specifically, a comparison with the reduc- 
tion factor of 3-D discs with Newtonian gravity reveals 
how well softening mimics the effects of thickness, as 
far as density waves are concerned. 

• At large \k\s, it shows how effectively noise is sup- 
pressed on small scales. 

• It also contains less direct information about how sig- 
nificantly the equilibrium properties are affected. 



The reduction factor 5(|fc|s) is related to the point-mass 
potential —Gm<p s (R) and force —Gm 2 f s (R) through inte- 
gral transforms. These transforms can easily be inverted, 
and the inversion formulae that relate <p s (R) and f s (R) to 
5(|fe|s) are: 



MR) = ^o 



S(\k\s) 



(R), 



f a {R)=H 1 [S(\k\a)](R), 

where TL U denotes the Hankcl transform of order v. 



n v [ g (k)\{R) 



g(k)J u (kR)kdk, 



(1) 



(2) 



(3) 



and J„ denotes the Bessel function of the first kind and or- 
der v (for mathematical and numerical references see Sect. 
II-2.1.1; in particular, for a swift numerical computation 
of Hankel transforms use, e.g., the NAG library). Inverting 
this part of the method has a direct practical importance 
and, indeed, is the trick behind the present application: it 
allows imposing arbitrary dynamical requirements in the 
form of a reduction factor and identifying the optimal type 
of softened gravity that satisfies such requirements. In ad- 
dition, the rest of the method allows testing the precise 
dynamical performance of the modelling. 

Here is an example: we want to find a type of soft- 
ened gravity that mimics the effects of thickness very well 
and, at the same time, suppresses noise very effectively. 
This is illustrated in Fig. la and in the following discus- 
sion. A reduction factor that conforms with the dynamical 
requirements specified above is: 



S(\k\s) 



l + \k\s 



-(l*l«/*) B 



(4) 



The first function is the reduction factor of 3-D discs with 
Newtonian gravity and characteristic scale height s (Shu 
1968; Vandervoort 1970; Romeo 1992). The second func- 
tion acts as a low-pass filter with rather sharp cut-off at 
radial wavelengths A « 2s. The benefits of filtering in spec- 
tral domain are well known in the context of digital image 
processing (see, e.g., Jain 1989; Press et al. 1992; see also 
Vetterling et al. 1992). In our case the dynamical reso- 
lution, i.e. the faithfulness in simulating the Newtonian 
dynamics, is substantially higher than in the standard 
Plummcr softening. Indeed, it can be further improved 
by choosing sharper filters, at the cost of noticeable os- 
cillations in the gravitational interaction. To some extent, 
their preference is a matter of taste. On the other hand, 
too sharp filters make the typical radial wavelength and 
possibly other dynamical properties hypersensitive to the 
location of the cut-off, which is unphysical. Thus, they 
should not be used. 
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2.2. Implementing star-gas softening 

How do wc model gravity in the presence of two com- 
ponents, such as stars and cold interstellar gas? Let us 
think in the alternative way of finite-sized particles inter- 
acting with Newtonian gravity, as Dyer & Ip (1993) have 
partly suggested (for an ABC of finite-sized particles see 
Appendix A). Then the answer is simple. Each component 
turns out to have its own positive reduction factor <S(|fc|s), 
where now s is the scale length of the particle mass dis- 
tribution. So the star-star and gas-gas interactions are as 
in the one-component case, while the star-gas interaction 
potential —Gm s m g ip s . g (R) and force —Gm s m g f s . g (R) are 
determined unequivocally by the inversion formulae: 



u s (R) = m 



- yfs B (\k\s B ) -S g (\k\s g ) 



(R). 



yJs B (\k\s B ) -S g (\k\s g 



(R). 



(5) 



(6) 



Our recipe for implementing star-gas softening has strong 
advantages. Indeed, its characteristics are fundamental for 



modelling the complex roles that such components play in 
regimes of astrophysical interest, as we have concluded in 
Paper I. 

As an example, we want to generalize the type of soft- 
ened gravity found in the one-component case to the pres- 
ence of a young disc stellar population and a cold inter- 
stellar gaseous component with, say, s s : s g — 2 : 1 (see, 
e.g., Mihalas & Binncy 1981). This is illustrated in Fig. 
lb and in the following discussion. The finite-sized particle 
implementation of star-gas softening is consistent with the 
effects of thickness: there are two positive reduction fac- 
tors, one for each component. Again, the stellar reduction 
factor is: 



5 s (|fc|s B ) = 



1 



-(|k|».A) B 



(7) 



Concerning cold interstellar gas, the situation is more 
complex and uncertain (e.g., Combes & Pfenniger 1996; 
Elmcgreen 1996a, b; Ferrara 1996; Lequeux & Guelin 
1996; Pfenniger 1996; Pfenniger et al. 1996). Sharp filters 
should not be used because they over-stress adherence to 
the effects of thickness, whereas the effects of turbulence 
and fractality may be more important. Soft filters are safer 
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Fig. 2. Examples of 3-D modelling: a 
isotropic case (cf. Sect. 3.1), b aniso- 
tropic case (cf. Sect. 3.2). The abbrevia- 
tions N, T and P mean Newtonian grav- 
ity, thickness and Plummer softening, re- 
spectively 



in that respect, and our preference goes to the Gaussian 
member of the familya. So the gaseous reduction factor is: 



" s ™-l + |fc| Sg 

Regarding s s as the characteristic scale height of refer- 
ence, different values of s g /s s have no influence on the 
stellar functions, modify the star-gas gravitational inter- 
action moderately and change the gaseous functions ac- 
cording to simple scaling laws. 

3. 3-D modelling 

3.1. Passing from 2D to 3D 

In 3-D disc models, <S(|fc|s) is no longer the true reduction 
factor but is still useful for quantifying the dynamical ef- 
fects of softening parallel to the plane, which now com- 
bine with those of vertical random motion. Specifically, 

2 A simpler, but less instructive, choice would be: <S g (|fc|s g ) = 
e~' fe ' Sg , i.e. the reduction factor for the standard Plummer 
softening. 



a comparison with the reduction factor of 3-D discs with 
Newtonian gravity suggests how much softening interferes 
with the effects of thickness, as far as density waves are 
concerned. The inversion formulae for <p s (r) and f s (r) are 
as in the 2-D case. This is the simplest way of passing 
from 2-D to 3-D modelling. (The generalization to two 
components is clear.) 

Here is an example: we want to find a type of softened 
gravity that interferes with the effects of thickness very 
little and, at the same time, suppresses noise very effec- 
tively. This is illustrated in Fig. 2a. The pseudo reduction 
factor is: 

5(|fc|s) = e- (|fc|s/ " )3 , (9) 

and acts as a low-pass filter with rather soft cut-off at 
A ~ 2s. The discussion concerning the dynamical resolu- 
tion and the action of sharper filters follows the 2-D case 
closely. On the other hand, only in 3D can we simulate 
the evolutionary nature of thickness, which arises from the 
vertical random motion and its subtle coupling with the 
dynamical properties parallel to the plane (Romeo 1990, 
1992; see also Paper I). 
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3.2. Implementing anisotropic softening 

In order to model gravity in situations that demand 
anisotropic resolution, it is convenient to think in the stan- 
dard way of point particles interacting with softened grav- 
ity, as Zotov & Morozov (1987) have partly suggested. The 
softening surface is transformed from a sphere of radius s 
into a spheroid of planar and vertical semi-axes sy and s±, 
respectively. This means that the softening length is the 
distance from the centre to the surface of the spheroid in 
the direction of the position vector: 

ls 2 R 2 + s 2 z 2 

«(frM) = V R 2 +Z 2 ■ (10) 

The resolution turns out to be decoupled parallel and per- 
pendicular to the plane, and to be determined by the natu- 
ral planar and vertical softening lengths. These character- 
istics of our recipe for implementing anisotropic softening 
have important advantages, as we have concluded in Paper 
II. (The generalization to two components is clear.) 

As an example, we want to generalize the type of 
softened gravity found in the isotropic case to situations 
that demand moderate anisotropic resolution with, say, 
S|| : s± = 2 : 1. This is illustrated in Fig. 2b. Regard- 
ing S|| as the softening length of reference, different values 
of s_i_/s|| have no influence on the planar functions and 
change the gravitational interaction along the vertical di- 
rection according to simple scaling laws. 

4. Discussion 

Three recent papers concerning softening optimization 
and conception deserve comment: 

• Merritt's (1996) optimization is performed with re- 
spect to the Newtonian dynamics in the configuration 
space, and concerns the softening length. The config- 
uration space does not permit a clear distinction be- 
tween large-scale dynamical properties and small-scale 
noise, and also emphasizes the equilibrium state as 
most representative of the whole dynamics. 

• Weinberg's (1996) optimization is comparable to that 
of Merritt (1996) but concerns orthogonal series force 
computation, i.e. roughly speaking the softening length 
and the type of softened gravity. 

• Dyer & Ip's (1993) conception is rigid: softened gravity 
must mimic finite-sized particles. But why? Softening 
is an artifice: its physical consistency should be scru- 
tinized with respect to basic dynamical requirements, 
not with respect to the inter-particle force alone. An 
elastic conception is more useful. Apart from that, 
Dyer & Ip (1993) have suggested a softening optimiza- 
tion that is not so different from that of Merritt (1996). 



5. Conclusions and perspectives 

The importance of computer simulations in astrophysics 
is analogous to that of experiments in other branches of 
physics. They also serve as a welcome bridge between the- 
ories, often restricted to idealized situations, and observa- 
tions, revealing instead the complexity of nature. Major 
present objectives are to construct physically consistent 
Af-body models of disc galaxies and to simulate their dy- 
namical evolution, especially in regimes of spiral struc- 
ture in which a fruitful comparison between theories and 
simulations can be made (e.g., Pfenniger & Friedli 1991; 
Junqueira & Combes 1996; Zhang 1996; Bottema & Ger- 
ritsen 1997; Fuchs & von Linden 1998; von Linden et al. 
1998; Zhang 1998a, b). The construction of such models 
is indeed a difficult task which has not yet been fully ac- 
complished, and which should eventually provide clues of 
vital importance to a number of open questions posed by 
both theories and observations. 

Our involvement has been threefold. In Paper I, we 
have recognized a fundamental problem posed by this re- 
search programme (for a concrete use of that analysis 
and for interesting remarks see, e.g., Junqueira & Combes 
1996). In Paper II, we have devised a method for solving 
this problem. In the present paper, we apply this method 
and solve the problem, thus laying the foundations of such 
a plan. The major result is that gravity can be modelled 
so as to optimize the fidelity of simulations, and the pro- 
cedure is practicable. The following conclusions point up 
the whys and wherefores: 

1. Optimization is performed with respect to arbitrary 
dynamical requirements and, in specific examples, with 
respect to the Newtonian dynamics. This enriches the 
modelling with an unprecedented degree of freedom, 
which has clear epistemological motivations (cf. Sect. 
1, discussion of the present paper). 

2. Optimization is performed in the wavenumber space. 
This is the appropriate domain for imposing dynamical 
requirements on the modelling. 

3. Optimization concerns both the softening length and 
the type of softened gravity. 

4. Softening is conceived as a double artifice. The soft- 
ened gravity and finite-sized particle conceptions are 
equivalent in the basic cases. Concerning more com- 
plex situations, the latter is particularly useful for im- 
plementing star-gas softening, whereas the former is 
particularly useful for implementing anisotropic soft- 
ening. Thus both conceptions contribute towards the 
accomplishment of a physically consistent modelling. 

Our application is ready for a concrete use. An attractive 
idea is to employ a particle-particle code together with 
MD-GRAPE, a highly parallelized special-purpose com- 
puter for many-body simulations with an arbitrary cen- 
tral force (Fukushige et al. 1996). We can also employ a 
classical particle-mesh code. Then the dynamical effects 
of the grid are known and factorize as those of soften- 
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ing (e.g., Bouchet et al. 1985; Efstathiou et al. 1985; for 
a review see, e.g., Hockney & Eastwood 1988). So essen- 
tially the application proceeds as in the present paper, 
but it may be useful to act directly on the wavenumber 
space (e.g., Tormen & Bertschinger 1996). A more com- 
plex problem concerns tree codes, which have hierarchi- 
cal structure and adaptive resolution over multiple scales 
(e.g., Hcrnquist 1987; for a review see, e.g., Pfalzner & 
Gibbon 1996). The solution to that problem would need 
a more advanced analysis (cf. following discussion). Wel- 
come suggestions about the choice of the code can come 
from cosmological simulations (e.g., Splinter et al. 1998). 

Finally, what about the future? Our approach is con- 
nected with the technique of filtering in spectral domain 
used in the context of digital image processing. This is 
a rapidly evolving field with growing applications in sci- 
ence and engineering, which can promote further sub- 
stantial advances in iV-body modelling of disc galaxies. 
For instance, wavelets are ideal for resolving multi-scale 
problems in space and/or time, such as those concerning 
turbulence, bifurcations, fractals and many others (see, 
e.g., Kaiser 1994; Holschncider 1995; Bowman & Newell 
1998; for an alternative analysis tool see, e.g., Stutzki et 
al. 1998). Speculating further, wavelets might be used for 
speeding up simulations through fast solution of linear sys- 
tems (cf. Press et al. 1992, pp. 597-599 and 782). 

These are the merits of our contribution. We hope 
that the trilogy (Papers I— III) and further reflections 
(Romeo 1998) will encourage A^-body experimenters to 
model gravity so as to optimize the fidelity of their sim- 
ulations, and that the result will be a stronger interdisci- 
plinary connection with theories and observations. 
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A. ABC of finite-sized particles 

Finite-sized particles interacting with Newtonian gravity 
are analogous to point particles interacting with softened 
gravity. The dynamics of one-component 2-D discs con- 
taining such particles can be investigated by performing 
an analysis comparable to that of Sects. II-2.1 and 2.1. 
In this appendix we report the formulae useful for Sect. 
2.2. Let m^, s (R) be the particle mass distribution of scale 
length s. The reduction factor is: 



The inversion formula for /i s (i?) is: 

f i s (R) = ^-H \VSW)} (R)- 
Zn L J 



(A2) 



Last and most useful, the inversion formulae for the inter- 
action potential —Gm 2 ^p s {R) and force —Gm 2 f s {R) are: 



ip s (R) = H 



■S(\k\s) 



f s (R)=n 1 [S(\k\s)](R) 



(A3) 



(A4) 



(For general references about finite-sized particles see, e.g., 
Vlasov 1961; Rohrlich 1965; Dawson 1983; Hockney & 
Eastwood 1988; Birdsall & Langdon 1991; in particular, 
the last reference contains useful insights in the context of 
A^-body simulations of plasmas.) 
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